<!DOCTYPE html>
<html>
<head><meta charset="utf-8" />
<title>Poissonfemrate</title><script src="https://cdnjs.cloudflare.com/ajax/libs/require.js/2.1.10/require.min.js"></script>
<script src="https://cdnjs.cloudflare.com/ajax/libs/jquery/2.0.3/jquery.min.js"></script>

<style type="text/css">
/* Overrides of notebook CSS for static HTML export */
body {
  overflow: visible;
  padding: 8px;
}
div#notebook {
  overflow: visible;
  border-top: none;
}@media print {
  div.cell {
    display: block;
    page-break-inside: avoid;
  } 
  div.output_wrapper { 
    display: block;
    page-break-inside: avoid; 
  }
  div.output { 
    display: block;
    page-break-inside: avoid; 
  }
}
</style>

<!-- Custom stylesheet, it must be in the parent directory as the html file -->
<link rel="stylesheet" type="text/css" media="all" href="../doc.css" />
<link rel="stylesheet" type="text/css" media="all" href="doc.css" />

<!-- Loading mathjax macro -->
<!-- Load mathjax -->
    <script src="https://cdnjs.cloudflare.com/ajax/libs/mathjax/2.7.1/MathJax.js?config=TeX-AMS_HTML"></script>
    <!-- MathJax configuration -->
    <script type="text/x-mathjax-config">
    MathJax.Hub.Config({
        tex2jax: {
            inlineMath: [ ['$','$'], ["\\(","\\)"] ],
            displayMath: [ ['$$','$$'], ["\\[","\\]"] ],
            processEscapes: true,
            processEnvironments: true
        },
        // Center justify equations in code and markdown cells. Elsewhere
        // we use CSS to left justify single line equations in code cells.
        displayAlign: 'center',
        "HTML-CSS": {
            styles: {'.MathJax_Display': {"margin": 0}},
            linebreaks: { automatic: true }
        }
    });
    </script>
    <!-- End of mathjax configuration --></head>
<body>
  <div tabindex="-1" id="notebook" class="border-box-sizing">
    <div class="container" id="notebook-container">

<div class="cell border-box-sizing text_cell rendered"><div class="prompt input_prompt">
</div>
<div class="inner_cell">
<div class="text_cell_render border-box-sizing rendered_html">
<h1 id="Linear-Element-for-Poisson-Equation-in-2D">Linear Element for Poisson Equation in 2D<a class="anchor-link" href="#Linear-Element-for-Poisson-Equation-in-2D">&#182;</a></h1>
</div>
</div>
</div>
<div class="cell border-box-sizing text_cell rendered"><div class="prompt input_prompt">
</div>
<div class="inner_cell">
<div class="text_cell_render border-box-sizing rendered_html">
<p>This example is to show the rate of convergence of the linear finite element approximation of the Poisson equation on the unit square:</p>
$$- \Delta u = f \; \hbox{in } (0,1)^2$$<p>for the following boundary conditions</p>
<ul>
<li>Non-empty Dirichlet boundary condition: $u=g_D \hbox{ on }\Gamma_D, \nabla u\cdot n=g_N \hbox{ on }\Gamma_N.$</li>
<li>Pure Neumann boundary condition: $\nabla u\cdot n=g_N \hbox{ on } \partial \Omega$.</li>
<li>Robin boundary condition: $g_R u + \nabla u\cdot n=g_N \hbox{ on }\partial \Omega$.</li>
</ul>
<p><strong>References</strong>:</p>
<ul>
<li><a href="femdoc.html">Quick Introduction to Finite Element Methods</a></li>
<li><a href="http://www.math.uci.edu/~chenlong/226/Ch2FEM.pdf">Introduction to Finite Element Methods</a></li>
<li><a href="http://www.math.uci.edu/~chenlong/226/Ch3FEMCode.pdf">Progamming of Finite Element Methods</a></li>
</ul>
<p><strong>Subroutines</strong>:</p>

<pre><code>- Poisson
- squarePoisson
- femPoisson
- Poissonfemrate

</code></pre>
<p>The method is implemented in <code>Poisson</code> subroutine and tested in <code>squarePoisson</code>. Together with other elements (P1, P2, P3, Q1), <code>femPoisson</code> provides a concise interface to solve Poisson equation. The P1 element is tested in <code>Poissonfemrate</code>. This doc is based on <code>Poissonfemrate</code>.</p>

</div>
</div>
</div>
<div class="cell border-box-sizing text_cell rendered"><div class="prompt input_prompt">
</div>
<div class="inner_cell">
<div class="text_cell_render border-box-sizing rendered_html">
<h2 id="P1-Linear-Element">P1 Linear Element<a class="anchor-link" href="#P1-Linear-Element">&#182;</a></h2><p>For the linear element on a simplex, the local basis functions are
barycentric coordinate of vertices. The local to global pointer is
<code>elem</code>. This is the simplest and default element for elliptic equations.</p>
<p><strong>A local basis of P1</strong></p>
<p>For $i = 1, 2, 3$, a local basis of the linear element space is given by the barycentric coordinate</p>
$$\phi_i = \lambda_i, \quad \nabla \phi_i = \nabla \lambda_i = - \frac{|e_i|}{2|T|}\boldsymbol n_i, \quad \int_T \nabla \phi_i\nabla \phi_j = \frac{1}{4|T|}\boldsymbol l_i \cdot \boldsymbol l_j$$<p>where $e_i$ is the edge opposite to the i-th vertex and $\boldsymbol n_i$ is the unit
outwards normal direction, and $\boldsymbol l_i$ is the edge vector of $e_i$.</p>
<p>See <a href="http://www.math.uci.edu/~chenlong/226/Ch2FEM.pdf">Finite Element Methods</a> Section 2.1 for geometric explanation of the barycentric coordinate and <a href="http://www.math.uci.edu/~chenlong/226/Ch3FEMcode.pdf">Programming of Finite Element Methods in MATLAB</a> for detailed explanation. For P1 element, the basic data structure <code>node,elem</code> is sufficient and displayed for the following coarse mesh.</p>

</div>
</div>
</div>
<div class="cell border-box-sizing code_cell rendered">
<div class="input">
<div class="prompt input_prompt">In&nbsp;[13]:</div>
<div class="inner_cell">
    <div class="input_area">
<div class=" highlight hl-matlab"><pre><span></span><span class="n">imatlab_export_fig</span><span class="p">(</span><span class="s">&#39;print-jpeg&#39;</span><span class="p">)</span>
<span class="n">node</span> <span class="p">=</span> <span class="p">[</span><span class="mi">0</span><span class="p">,</span><span class="mi">0</span><span class="p">;</span> <span class="mi">1</span><span class="p">,</span><span class="mi">0</span><span class="p">;</span> <span class="mi">1</span><span class="p">,</span><span class="mi">1</span><span class="p">;</span> <span class="mi">0</span><span class="p">,</span><span class="mi">1</span><span class="p">];</span>
<span class="n">elem</span> <span class="p">=</span> <span class="p">[</span><span class="mi">2</span><span class="p">,</span><span class="mi">3</span><span class="p">,</span><span class="mi">1</span><span class="p">;</span> <span class="mi">4</span><span class="p">,</span><span class="mi">1</span><span class="p">,</span><span class="mi">3</span><span class="p">];</span>      
<span class="p">[</span><span class="n">node</span><span class="p">,</span><span class="n">elem</span><span class="p">]</span> <span class="p">=</span> <span class="n">uniformbisect</span><span class="p">(</span><span class="n">node</span><span class="p">,</span><span class="n">elem</span><span class="p">);</span>
<span class="n">figure</span><span class="p">(</span><span class="s">&#39;rend&#39;</span><span class="p">,</span><span class="s">&#39;painters&#39;</span><span class="p">,</span><span class="s">&#39;pos&#39;</span><span class="p">,[</span><span class="mi">10</span> <span class="mi">10</span> <span class="mi">225</span> <span class="mi">225</span><span class="p">])</span>
<span class="n">showmesh</span><span class="p">(</span><span class="n">node</span><span class="p">,</span><span class="n">elem</span><span class="p">);</span>
<span class="n">findnode</span><span class="p">(</span><span class="n">node</span><span class="p">);</span>
<span class="n">findelem</span><span class="p">(</span><span class="n">node</span><span class="p">,</span><span class="n">elem</span><span class="p">);</span>
<span class="n">display</span><span class="p">(</span><span class="n">elem</span><span class="p">);</span>
</pre></div>

</div>
</div>
</div>

<div class="output_wrapper">
<div class="output">


<div class="output_area">

<div class="prompt"></div>


<div class="output_subarea output_stream output_stdout output_text">
<pre>
elem =

     8     6     2
     7     6     4
     5     6     1
     9     6     3
     8     3     6
     7     1     6
     5     2     6
     9     4     6

</pre>
</div>
</div>

<div class="output_area">

<div class="prompt"></div>




<div class="output_jpeg output_subarea ">
<img src=""
>
</div>

</div>

</div>
</div>

</div>
<div class="cell border-box-sizing code_cell rendered">
<div class="input">
<div class="prompt input_prompt">In&nbsp;[14]:</div>
<div class="inner_cell">
    <div class="input_area">
<div class=" highlight hl-matlab"><pre><span></span><span class="c">% Setting</span>
<span class="p">[</span><span class="n">node</span><span class="p">,</span><span class="n">elem</span><span class="p">]</span> <span class="p">=</span> <span class="n">squaremesh</span><span class="p">([</span><span class="mi">0</span><span class="p">,</span><span class="mi">1</span><span class="p">,</span><span class="mi">0</span><span class="p">,</span><span class="mi">1</span><span class="p">],</span><span class="mf">0.25</span><span class="p">);</span> 
<span class="n">mesh</span> <span class="p">=</span> <span class="n">struct</span><span class="p">(</span><span class="s">&#39;node&#39;</span><span class="p">,</span><span class="n">node</span><span class="p">,</span><span class="s">&#39;elem&#39;</span><span class="p">,</span><span class="n">elem</span><span class="p">);</span>
<span class="n">figure</span><span class="p">(</span><span class="s">&#39;rend&#39;</span><span class="p">,</span><span class="s">&#39;painters&#39;</span><span class="p">,</span><span class="s">&#39;pos&#39;</span><span class="p">,[</span><span class="mi">10</span> <span class="mi">10</span> <span class="mi">225</span> <span class="mi">225</span><span class="p">])</span>
<span class="n">showmesh</span><span class="p">(</span><span class="n">node</span><span class="p">,</span><span class="n">elem</span><span class="p">);</span>
<span class="n">option</span><span class="p">.</span><span class="n">L0</span> <span class="p">=</span> <span class="mi">3</span><span class="p">;</span>
<span class="n">option</span><span class="p">.</span><span class="n">maxIt</span> <span class="p">=</span> <span class="mi">4</span><span class="p">;</span>
<span class="n">option</span><span class="p">.</span><span class="n">printlevel</span> <span class="p">=</span> <span class="mi">1</span><span class="p">;</span>
</pre></div>

</div>
</div>
</div>

<div class="output_wrapper">
<div class="output">


<div class="output_area">

<div class="prompt"></div>




<div class="output_jpeg output_subarea ">
<img src=""
>
</div>

</div>

</div>
</div>

</div>
<div class="cell border-box-sizing text_cell rendered"><div class="prompt input_prompt">
</div>
<div class="inner_cell">
<div class="text_cell_render border-box-sizing rendered_html">
<h2 id="Mixed-boundary-condition">Mixed boundary condition<a class="anchor-link" href="#Mixed-boundary-condition">&#182;</a></h2>
</div>
</div>
</div>
<div class="cell border-box-sizing code_cell rendered">
<div class="input">
<div class="prompt input_prompt">In&nbsp;[12]:</div>
<div class="inner_cell">
    <div class="input_area">
<div class=" highlight hl-matlab"><pre><span></span><span class="n">option</span><span class="p">.</span><span class="n">plotflag</span> <span class="p">=</span> <span class="mi">0</span><span class="p">;</span>
<span class="n">pde</span> <span class="p">=</span> <span class="n">sincosdata</span><span class="p">;</span>
<span class="n">mesh</span><span class="p">.</span><span class="n">bdFlag</span> <span class="p">=</span> <span class="n">setboundary</span><span class="p">(</span><span class="n">node</span><span class="p">,</span><span class="n">elem</span><span class="p">,</span><span class="s">&#39;Dirichlet&#39;</span><span class="p">,</span><span class="s">&#39;~(x==0)&#39;</span><span class="p">,</span><span class="s">&#39;Neumann&#39;</span><span class="p">,</span><span class="s">&#39;x==0&#39;</span><span class="p">);</span>
<span class="n">femPoisson</span><span class="p">(</span><span class="n">mesh</span><span class="p">,</span><span class="n">pde</span><span class="p">,</span><span class="n">option</span><span class="p">);</span>
</pre></div>

</div>
</div>
</div>

<div class="output_wrapper">
<div class="output">


<div class="output_area">

<div class="prompt"></div>


<div class="output_subarea output_stream output_stdout output_text">
<pre>Warning: File: /Dropbox/Math/Programming/ifem/solver/mg.m Line: 695 Column: 24
Defining &#34;Di&#34; in the nested function shares it with the parent function.  In a future release, to share &#34;Di&#34; between parent and nested functions, explicitly define it in the parent function.
&gt; In Poisson (line 233)
  In femPoisson (line 65)
Multigrid V-cycle Preconditioner with Conjugate Gradient Method
#dof:     4225,  #nnz:    19906, smoothing: (1,1), iter: 10,   err = 1.55e-09,   time = 0.11 s
Multigrid V-cycle Preconditioner with Conjugate Gradient Method
#dof:    16641,  #nnz:    80770, smoothing: (1,1), iter: 10,   err = 1.55e-09,   time = 0.096 s
Multigrid V-cycle Preconditioner with Conjugate Gradient Method
#dof:    66049,  #nnz:   325378, smoothing: (1,1), iter: 10,   err = 1.47e-09,   time =  0.2 s

 #Dof       h       ||u-u_h||    ||Du-Du_h||   ||DuI-Du_h|| ||uI-u_h||_{max}

 1089   3.12e-02   1.15027e-03   1.08974e-01   2.21506e-03   9.04547e-04
 4225   1.56e-02   2.88013e-04   5.45135e-02   5.54571e-04   2.26928e-04
16641   7.81e-03   7.20310e-05   2.72601e-02   1.38693e-04   5.67600e-05
66049   3.91e-03   1.80095e-05   1.36305e-02   3.46767e-05   1.41918e-05

 #Dof   Assemble     Solve      Error      Mesh    

 1089   7.00e-02   1.36e-02   9.00e-02   1.00e-02
 4225   3.00e-02   1.09e-01   4.00e-02   2.00e-02
16641   1.10e-01   9.62e-02   8.00e-02   6.00e-02
66049   4.00e-01   2.00e-01   2.30e-01   2.40e-01


</pre>
</div>
</div>

<div class="output_area">

<div class="prompt"></div>




<div class="output_jpeg output_subarea ">
<img src=""
>
</div>

</div>

</div>
</div>

</div>
<div class="cell border-box-sizing text_cell rendered"><div class="prompt input_prompt">
</div>
<div class="inner_cell">
<div class="text_cell_render border-box-sizing rendered_html">
<h2 id="Pure-Neumann-boundary-condition">Pure Neumann boundary condition<a class="anchor-link" href="#Pure-Neumann-boundary-condition">&#182;</a></h2><p>When pure Neumann boundary condition is posed, i.e., $-\Delta u =f$ in $\Omega$ and $\nabla u\cdot n=g_N$ on $\partial \Omega$, the data should be consisitent in the sense that $\int_{\Omega} f \, dx + \int_{\partial \Omega} g \, ds = 0$. The solution is unique up to a constant. A post-process is applied such that the constraint $\int_{\Omega}u_h dx = 0$ is imposed.</p>

</div>
</div>
</div>
<div class="cell border-box-sizing code_cell rendered">
<div class="input">
<div class="prompt input_prompt">In&nbsp;[13]:</div>
<div class="inner_cell">
    <div class="input_area">
<div class=" highlight hl-matlab"><pre><span></span><span class="n">option</span><span class="p">.</span><span class="n">plotflag</span> <span class="p">=</span> <span class="mi">0</span><span class="p">;</span>
<span class="n">pde</span> <span class="p">=</span> <span class="n">sincosNeumanndata</span><span class="p">;</span>
<span class="n">pde</span> <span class="p">=</span> <span class="n">sincosdata</span><span class="p">;</span>
<span class="n">mesh</span><span class="p">.</span><span class="n">bdFlag</span> <span class="p">=</span> <span class="n">setboundary</span><span class="p">(</span><span class="n">node</span><span class="p">,</span><span class="n">elem</span><span class="p">,</span><span class="s">&#39;Neumann&#39;</span><span class="p">);</span>
<span class="n">femPoisson</span><span class="p">(</span><span class="n">mesh</span><span class="p">,</span><span class="n">pde</span><span class="p">,</span><span class="n">option</span><span class="p">);</span>
</pre></div>

</div>
</div>
</div>

<div class="output_wrapper">
<div class="output">


<div class="output_area">

<div class="prompt"></div>


<div class="output_subarea output_stream output_stdout output_text">
<pre>Multigrid V-cycle Preconditioner with Conjugate Gradient Method
#dof:     4225,  #nnz:    20860, smoothing: (1,1), iter: 12,   err = 1.43e-09,   time = 0.03 s
Multigrid V-cycle Preconditioner with Conjugate Gradient Method
#dof:    16641,  #nnz:    82684, smoothing: (1,1), iter: 12,   err = 6.33e-09,   time = 0.08 s
Multigrid V-cycle Preconditioner with Conjugate Gradient Method
#dof:    66049,  #nnz:   329212, smoothing: (1,1), iter: 13,   err = 2.40e-09,   time = 0.26 s

 #Dof       h       ||u-u_h||    ||Du-Du_h||   ||DuI-Du_h|| ||uI-u_h||_{max}

 1089   3.12e-02   1.29973e-03   1.08855e-01   5.80361e-03   3.86104e-03
 4225   1.56e-02   3.25931e-04   5.44960e-02   1.64867e-03   1.14414e-03
16641   7.81e-03   8.15520e-05   2.72576e-02   4.64515e-04   3.30465e-04
66049   3.91e-03   2.03927e-05   1.36301e-02   1.29765e-04   9.37017e-05

 #Dof   Assemble     Solve      Error      Mesh    

 1089   5.00e-02   1.49e-03   2.00e-02   1.00e-02
 4225   2.00e-02   3.03e-02   3.00e-02   3.00e-02
16641   8.00e-02   7.97e-02   6.00e-02   5.00e-02
66049   4.00e-01   2.65e-01   3.00e-01   2.60e-01


</pre>
</div>
</div>

<div class="output_area">

<div class="prompt"></div>




<div class="output_jpeg output_subarea ">
<img src=""
>
</div>

</div>

</div>
</div>

</div>
<div class="cell border-box-sizing text_cell rendered"><div class="prompt input_prompt">
</div>
<div class="inner_cell">
<div class="text_cell_render border-box-sizing rendered_html">
<h2 id="Robin-boundary-condition">Robin boundary condition<a class="anchor-link" href="#Robin-boundary-condition">&#182;</a></h2>
</div>
</div>
</div>
<div class="cell border-box-sizing code_cell rendered">
<div class="input">
<div class="prompt input_prompt">In&nbsp;[14]:</div>
<div class="inner_cell">
    <div class="input_area">
<div class=" highlight hl-matlab"><pre><span></span><span class="n">option</span><span class="p">.</span><span class="n">plotflag</span> <span class="p">=</span> <span class="mi">0</span><span class="p">;</span>
<span class="n">pde</span> <span class="p">=</span> <span class="n">sincosRobindata</span><span class="p">;</span>
<span class="n">mesh</span><span class="p">.</span><span class="n">bdFlag</span> <span class="p">=</span> <span class="n">setboundary</span><span class="p">(</span><span class="n">node</span><span class="p">,</span><span class="n">elem</span><span class="p">,</span><span class="s">&#39;Robin&#39;</span><span class="p">);</span>
<span class="n">femPoisson</span><span class="p">(</span><span class="n">mesh</span><span class="p">,</span><span class="n">pde</span><span class="p">,</span><span class="n">option</span><span class="p">);</span>
</pre></div>

</div>
</div>
</div>

<div class="output_wrapper">
<div class="output">


<div class="output_area">

<div class="prompt"></div>


<div class="output_subarea output_stream output_stdout output_text">
<pre>Multigrid V-cycle Preconditioner with Conjugate Gradient Method
#dof:     4225,  #nnz:    20865, smoothing: (1,1), iter: 10,   err = 1.99e-09,   time = 0.018 s
Multigrid V-cycle Preconditioner with Conjugate Gradient Method
#dof:    16641,  #nnz:    82689, smoothing: (1,1), iter: 10,   err = 2.33e-09,   time = 0.062 s
Multigrid V-cycle Preconditioner with Conjugate Gradient Method
#dof:    66049,  #nnz:   329217, smoothing: (1,1), iter: 10,   err = 3.17e-09,   time = 0.21 s

 #Dof       h       ||u-u_h||    ||Du-Du_h||   ||DuI-Du_h|| ||uI-u_h||_{max}

 1089   3.12e-02   4.92975e-03   4.34581e-01   2.56571e-02   8.30859e-03
 4225   1.56e-02   1.24034e-03   2.17889e-01   6.44198e-03   2.08620e-03
16641   7.81e-03   3.10581e-04   1.09020e-01   1.61223e-03   5.22032e-04
66049   3.91e-03   7.76764e-05   5.45192e-02   4.03168e-04   1.30532e-04

 #Dof   Assemble     Solve      Error      Mesh    

 1089   5.00e-02   1.39e-03   0.00e+00   1.00e-02
 4225   2.00e-02   1.76e-02   2.00e-02   3.00e-02
16641   8.00e-02   6.17e-02   6.00e-02   5.00e-02
66049   4.90e-01   2.13e-01   2.50e-01   2.40e-01


</pre>
</div>
</div>

<div class="output_area">

<div class="prompt"></div>




<div class="output_jpeg output_subarea ">
<img src=""
>
</div>

</div>

</div>
</div>

</div>
<div class="cell border-box-sizing text_cell rendered"><div class="prompt input_prompt">
</div>
<div class="inner_cell">
<div class="text_cell_render border-box-sizing rendered_html">
<h2 id="Conclusion">Conclusion<a class="anchor-link" href="#Conclusion">&#182;</a></h2><p>The optimal rate of convergence of the H1-norm (1st order) and L2-norm
(2nd order) is observed. The 2nd order convergent rate between two
discrete functions $\|\nabla (u_I - u_h)\|$ is known as superconvergence.</p>
<p>MGCG converges uniformly in all cases.</p>

</div>
</div>
</div>
    </div>
  </div>
</body>

 


</html>
